function [P, Cost ] = optimizeFeatureProgram( X, P, W, Y, MAXIT, alpha, LB, WB )
%optimizeFeatureProgram is a linear regression to optimize the feature
%program parameters
%  X is the raw data, Y is the target feature, MAXIT is the number of iterations, 
%Alpha is the learning rate and lambda is the regularization factor
%RETURN: Parameters P, Cost is the cost function for each iteration

    
    %number of features
    N = size(X,2);
    
    %transform the p parameters in a column vector
    P = P(:);
    
    %Store the previous P
    PreviousP = P(:);
    
    %create a vector to store the cost function in each iteration
    Cost = zeros(MAXIT,1);
        
    %loop for each iterarion
    for i = 1:MAXIT
    

        %calculate new features
        [FY, F] = featureProgram(Y, X, reshape(P,N,2));
        
        %number of examples
        M = size(F,1);
        
        %calculate the hypothesis. obs: intercept term
        H = [ones(M,1) F] * W; 
        
        %calculate the cost
        Cost(i,1) = (1/(2*M)) .* ( sum((H - FY) .^ 2));
        
        if i == 1
            %initial pertubation. 2n because p was nX2 matrix
            P = P + sign((-1 + (1+1) .* rand(2*N,1)));
            
            %constraint the new parameters to positive
            P =((P >= 0) .* P) + (P == 0);
            
        else
            %adjustment based in how much the erro change according to the
            %direction in P adjustment
            costDir = sign(Cost(i,1) - Cost((i - 1),1));
            if isnan(costDir)
                costDir = 1;
            end;
            
            %calculate the direction
            parDir = sign(P - PreviousP);
            
            %save the previous P to the next iteration
            PreviousP = P;
            
            %adjust the parameters
            P = round(P - (alpha .* parDir .* costDir));
            
            %constraint the new parameters to positive
            P =((P >= 0) .* P) + (P == 0);
        end;
        
    end;
    
    %reshape the p parameters to Mx2 matrix
    P = reshape(P,N,2);
end

